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We review the Color Glass Condensate effective theory, that describes the gluon content 
of a high energy hadron or nucleus, in the saturation regime. The emphasis is put on ap- 
plications to high energy heavy ion collisions. After describing initial state factorization, 
we discuss the Glasma phase, that precedes the formation of an equilibrated quark-gluon 
plasma. We end this review with a presentation of recent developments in the study of 
the isotropization and thermalization of the quark-gluon plasma. 

Keywords: Heavy Ion collisions, Quantum Chromodynamics, Gluon saturation, Color 
Glass Condensate, Glasma, Thermalization 

PACS numbers: 12.20.-ni, ll.10.Hi, 11.15.Kc, ll.80.Fv, 11. 80. La, 13.85.Hd 
1. Introduction 

Heavy Ion Collisions at high energy are performed at the Relativistic Heavy Ion Col- 
lider (RHIC) and the Large Hadron Collider (LHC), in order to study the properties 
of nuclear matter under extreme conditions of density and temperature. From lattice 
simulations,^ it is known that above a certain critical temperature, the protons and 
neutrons should melt into a plasma made of their constituents, the Quark-Gluon 
Plasma (QGP), and the critical temperature of this deconfinement transition is 
expected to be reached in the collisions performed at RHIC and LHC. 




Fig. 1. Elementary couplings between quarks and gluons in QCD. 

Quantum Chromodynamics (QCD), the theory of the microscopic interactions 

*Based on lectures given at the 22nd Jyvaskyla Summer School, August 2012, Jyvaskyla, Finland. 
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between the quarks and the gluons (see the figure [T]) , should in principle also be 
applicable to heavy ion collisions. However, since the strong interaction coupling 
constant becomes large at low momentum, it is not obvious a priori that heavy ion 
collisions can be studied by weak coupling techniques. This is certainly possible for 
the rare large-momentum processes that take place in these collisions, but quite 
questionable for the bulk of the particle production processes. Moreover, since the 
system produced in such a collision expands rapidly along the collision axis (see 
the figure [2]), its characteristic momentum scales (e.g. its temperature if it reaches 
thermal equilibrium) decrease with time. Therefore, there will always be a time 



Fig. 2. Successive stages of a high energy heavy ion collision, and the most widely used approaches 
to describe each stage. 

beyond which the coupling is strong and weak coupling approaches are useless. 
This is obviously the case near the phase transition that sees the quarks and gluons 
recombine in order to form hadrons. In the best of cases, we can only hope for a 
weak coupling treatment of the early stages of these collisions (say up to a couple 
of fm/c after the collision). 



Fig. 3. Parton content of a nucleon or nucleus at low energy. Left: cartoon of a nucleus at low 
energy and its valence quarks. Right: the thick lines represent the valence quarks, and the wavy 
lines are gluons. Virtual quark-antiquark pairs are not represented. 



t 





When applying QCD to the study of hadronic collisions, an essential ingredi- 
ent is the quark and gluon content of the hadrons that are being collided, since 
the elementary degrees of freedom in QCD arc partons rather than hadrons. On the 
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surface, a rmcleon is made of three valence quarks, bound by gluons. However, these 
quarks can also temporarily fluctuate into states that have additional gluons and 
quark-antiquark pairs (see the figure [3]). These fluctuations are short lived, with a 
lifetime that is inversely proportional to their energy. The largest possible lifetime 
of these fluctuations is comparable to the nucleon size, and they can be arbitrarily 
short lived. However, in a given reaction that probes the nucleon, there is always a 
characteristic time scale set by the resolution power of the probe (for instance by 
the frequency of the virtual photon that probes the nucleon in Deep Inelastic Scat- 
tering) . Only the fluctuations that are longer lived than the resolution in time of the 
probe can actually be seen in the process. The shorter lived fluctuations are present, 
but do not influence the reaction. In collisions involving a low energy nucleon, only 
its valence quarks and a few of these fluctuations are visible. Moreover, in a low 
energy nucleon, there will typically be interactions between its constituents during 
the collision with the probe, thus making low energy reactions very complicated. 




Fig. 4. Parton content of a nucleon or nucleus at high energy, as seen in the laboratory frame. 
Left: boosted nucleus and its gluonic content. Right: fluctuations inside a boosted nucleon. 

However, this picture is dramatically modified when the reaction involves a high 
energy nucleon, due to relativistic kinematics (see the figure[4|. Firstly, the geometry 
of the nucleon changes due to Lorentz contraction: at very high energy, the nucleon 
appears essentially two-dimensional in the laboratory fram^j] Simultaneously, all 
the internal timcscalcs of the nucleon in particular the lifetimes of the fluctuations 
and the duration of the interactions among the constituents- are multiplied by 
the same Lorentz factor. The first consequence of this time dilation is that the 
partons are now unlikely to interact precisely during the time interval probed in 
the reaction: the constituents of a high energy nucleon appear to be free during 
the collision. Secondly, since the lifetimes of the fluctuations are also dilated, more 
fluctuations are now visible by the probe: the number of gluons seen in a reaction 
increases with the energy of the collision. This increase with energy of the number of 
gluons in a nucleon has been observed experimentally in Deep Inelastic Scattering 
(DIS), for instance at HERA. This is shown in the figure [5] for a proton. Note that 
in this plot, high energy corresponds to small values of the longitudinal momentum 
fraction x carried by the parton, x = p z /y/s. The other important feature of the 

a The Lorentz gamma factor is 7 ~ 100 at RHIC and 7 ~ 1000 in heavy ion collisions at the LHC. 
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Fig. 5. Parton distributions in a proton, measured in Deep Inelastic Scattering at HERA. 2 



parton distributions at high energy/small x is that the gluons are outnumbering 
all the other parton species - the valence quarks are completely negligible in this 
kinematical region, and the sea quarks are suppressed by one power of the coupling 
a s , since they are produced from the gluons by the splitting g — > qq. 




Fig. 6. Left: typical process in a hadronic collision in the dilute regime. Right: typical process in 
the dense regime. 



This increase of the gluon distribution at small x leads to a major complication 
when applying QCD to compute processes in this regime. Indeed, the usual tools of 
perturbation theory are well adapted to the situation where the parton distributions 
are small (see the left figure [6| and where a fairly small number of graphs contribute 
at each order. On the contrary, when the parton distributions increase, processes 
involving many partons become more and more important, as illustrated in the 
right panel of the figure [6] The extreme situation arises when the gluon occupation 
number is of order l/a s : in this case, an infinite number of graphs contribute at each 
order. This regime of high parton densities is non-perturbative, even if the coupling 
constant is weak - the non-perturbative features arise from the fact that the large 
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parton density compensates the smallness of the coupling constant. In addition 
to requiring the summation of an infinite set of Feynman diagrams, this situation 
requires some knowledge about the probability of occurrence of multigluon states in 
the wavefunctions of the two colliding projectiles, something which is not provided 
by the usual parton distributions (they give only the single parton density). 

A major progress in dealing with this situation has been the realization that weak 
coupling methods can be used in this problem, thanks to the dynamical generation 
of a scale that is much larger than the non-perturbative scale A QCD (the scale at 
which the strong interactions become really strong, and where confinement effects 
are crucial). This scale, known as the saturation momentum and denoted Q s , is due 
to the non-linear interactions among the gluons. Roughly speaking, it is defined as 
the coupling constant a s times the gluon density per unit area (because the Lorentz 
contraction makes the nucleus look like a flat sheet of gluons in the laboratory 
frame), and it is a measure of the strength of the gluon recombination processes 
that may occur when the gluon density becomes largeP^ Any process involving 
momenta smaller than Q s may be affected by gluon saturation. 

Moreover, an effective theory -known as the Color Glass Condensate (CGC)- has 
been developed in order to organize the calculations of processes in the saturation 
regime. This effective theory approximates the description of the fast partons in the 
wavefunction of a hadron by exploiting the fact that their dynamics is slowed down 
by Lorentz time dilation, and provides a way to track the evolution with energy 
of the multigluon states that are relevant in the dense regime. This framework has 
been applied to a range of reactions at high energy: DIS, proton-nucleus collisions 
and nucleus-nucleus collisions. At leading order, these calculations correspond to a 
classical field description of the system. In nucleus-nucleus collisions, this classical 
field remains coherent for a brief amount of time after the collision, forming a state 
know as the Glasma. A central question in heavy ion collisions is to understand how 
these classical fields lose their coherence in order to form a plasma of quarks and 
gluons in local thermal equilibrium. 

The purpose of these lectures is to expose the physics of the Color Glass Con- 
densate and of the Glasma. The main focus will be heavy ion collision^] except for 
the first sections where the simpler example of DIS is used in order to introduce the 
concept of gluon saturation and the CGC framework. Our emphasis is to present a 
consistent framework to describe heavy ion collisions from the pre-collision initial 
state to the time at which the system may be described by hydrodynamics]^] The 
outline is the following: 

• In the section [2j we discuss the evolution of the wavefunction of a color 
singlet dipole, and derive the BFKL and BK equations. This will provide a 

b For this reason, the topic of pomeron loops will not be covered here, because these effects are 
important only in the dilute regimePHUl 

c Theoretical and phenomcnological aspects of the Color Glass Condensate, with slightly different 
points of view, have also been addressed in other reviewsPSGH 
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first glimpse at gluon saturation, and allow us to explain some interesting 
scaling properties for the inclusive DIS cross-section. 

• In the section [3j we introduce the Color Glass Condensate effective theory, 
as a more general way to describe the saturation regime. 

• The section [4] is devoted to a discussion of the factorization of the loga- 
rithms that appear in loop calculations in the CGC framework, for heavy 
ion collisions. We show how these logarithms can be absorbed into univer- 
sal distributions that describe the multigluon content of the two colliding 
nuclei. 

• In the section [5j we present the main ideas of the Glasma picture, and 
discuss some phenomenological consequences that follow from it. 

• The section [6] discusses the final state evolution, that leads from the coher- 
ent Glasma fields to the thermalized quark gluon plasma. The main focus 
is on the instabilities of the classical solutions of the Yang-Mills equations, 
their resummation, and their possible role in the isotropization and ther- 
malization processes. 

2. Hadron wavefunction at high energy 

In order to introduce the physics of gluon saturation, we consider in this section 
the case of Deep Inelastic Scattering. This situation is one of the simplest, since 
only one of the two projectiles is a hadron, the other being a lepton that interacts 
only electromagnetically via a photon exchange. DIS can viewed as an interaction 
between a photon with a negative virtuality {q^q^ = —Q 2 < 0) and the hadronic 
target. 




Fig. 7. Cartoon of Deep Inelastic Scattering in a frame where the photon fluctuates into a dipole. 

The cross-section for this process is frame independent, but in this section we 
will analyze it in a frame where the photon has a large longitudinal momentum, 
while the hadron has only a moderate momentum. In such a frame, the photon 
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can fluctuate into a quark-antiquark paiij^] in a color singlet state (called a color 
dipole), and it is this pair that interacts with the colored constituents of the target. 
Moreover, at high energy, this dipole sees the constituents (quarks and gluons) 
as frozen due to time dilation. Therefore, its interactions with these constituents 
can be approximated as interactions with the (static) color field that they create. 
Moreover, thanks to confinement, it is legitimate to assume that this target color 
field occupies a bounded region of space. 



2.1. Eikonal approximation 

We will compute the total cross-section between the dipole and the target field by 
using the optical theorem, that relates the total cross-section to the forward scat- 
tering amplitude. In the limit where the longitudinal momentum of the dipole is 
very large, this amplitude becomes very simple and is given by the eikonal approx- 
imation. Since this is an extremely important result in the study of high energy 
scattering, let us derive it in detail. Consider an S'-matrix element, 

S/3a = (0oixt\oLva) = (An| U(+OC, -Oo) \of m ) , (1) 

for the transition between two arbitrary states made of quarks, antiquarks and 
gluons, a and /?. In the second equality, C/(+oo, — oo) is the evolution operator from 
the initial to the final state. It can be expressed as the time ordered exponential of 
the interaction part of the Lagrangian, 

U(+oo, — oo) = T exp i J d 4 x £ int ((j)i n (x)) , (2) 

where <p- m denotes generically the fields in the interaction picture. In our problem, 
Ant contains both the self-interactions of the fields, and their interaction with the 
target color field. We want to compute the high energy limit of this scattering 
amplitude, 

sfc> = Jin^ <An|e-^V(+oo, -oo) e+^ 3 |a in ) (3) 

boosted state 

where K 3 is the generator of Lorentz boosts in the +z direction. 

Before doing any calculation, a simple argument can help understand what hap- 
pens in this limit. Quite generally, scattering amplitudes arc proportional to the 
duration of the overlap between the wavefunctions of the two colliding objects. In 
the present case, it should scale as the time spent by the incoming state in region 
occupied by the target field. This time is inversely proportional to the energy of the 
incoming state, and goes to zero in the limit uj — > +oo. If the interaction between the 
projectile and the target field was via a scalar exchange, then the conclusion would 



d Of course, the photon may also fluctuate into more complicated Fock states, such as a qqg state, 
but the probability of occurrence of these states is suppressed by at least one extra power of the 
coupling a s and therefore they do not contribute at leading order. 
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be that the scattering amplitude vanishes in the high energy limit (in other words, 
S'-matrix elements would go to unity). However, interactions with a color field are 
via a vector exchange, i.e. the target field couples to a four-vector that represents 
the color current carried by the projectile, by a term of the form —igA^J^. At high 
energy, the longitudinal component of this four-vector increases proportionally to 
the energy, and compensates the small time spent in the interaction zone. Thus, for 
states that interact via a vector exchange^ we expect that scattering amplitudes 
have a finite high energy limit (nor zero, nor infinite). 

This calculation is best done using light-cone coordinates. For any four-vector 
a p , one defines 



V2 ' ~ V2 
The following formulas are often useful, 



(4) 



x ■ y = x + y + x y + - x± y ± 
d A x = dx + dx~ d 2 x±^ 

□ = 2d+d- - Vi with 8+ = — — , d~ = t— . (5) 

ox~ ox + 

For a highly boosted projectile in the +z direction, x + plays the role of the time, 
and the Hamiltonian is the P~ component of the momentum. The generator of 
longitudinal boosts in light-cone coordinates is 

R- = M~ + = -K 3 . (6) 

In order to derive the eikonal approximation, the following identities are also very 
useful, 

e iuK~ p+ ^ — iu:K~ _ e + u > p+ 

e iwK- pj e -iuK~ = pj _ ^) 

They express the fact that, under longitudinal boosts, the components P ± of a four- 
vector are simply rescaled, while the transverse components are left unchanged. 

Under longitudinal boosts, states, creation operators and field operators are 
transformed as follows, 

e- iuK ~\p---i a ) = \(e"p+, Px ) ■■■ iB ) 

e~ iuK ' al(q)e^ K ~ = 4 a (e"q + ,e-»q-,q ± ) 

e iuK ~ <f> in (x)e- iuK ~ = <l> in (e- u x+,e" 3 x-,x ± ) . (8) 

Note that the last equation is valid only for a scalar field, or for the transverse 
components of a vector field. In addition, the ± components of a vector field receive 



c By the same reasoning, gravitational interactions, that involve a spin two exchange, would lead 
to scattering amplitudes that grow linearly with energy. 



Color Glass Condensate and Glasma 9 



an overall rescaling by a factor e ±LO . Moreover, since a longitudinal boost does not 
alter the time ordering, we can also write 



?7(+oo, — oo)e 



Texpi I d*x £ int (e luK ~ <p in (x)e 



)• (9) 



Likewise, the components of the vector current that couples to the target field 
transform as 



j\x)e 



iuK J-{x)e- iuK ~ 



= J i {e-° J x + ,e 0J x-,x ± ) 
= e" w J-{e~ w x + ,e u x- ,xj_) 



J + (x)e 



" - e u J+ie'^x+^e^x-^x^) . 



(10) 



Naturally, the target field A^ does not change when we boost the projectile. For 
simplicity, let us assume that A^ is confined in the region — L < x + < +L. We can 
split the evolution operator into three factors, 



U(+oo, -oo) = U(+oo, +L)U(+L, -L)U{-L, -oo) 



(11) 



The factors f/(+oo, +L) and U(—L, — oo) do not contain the external potential. For 
these two factors, the change of variables e~ u x + — > x + , e u x~ — > x~ leads to 



lim e wK U(+oo,+L)e 

oj—>-\-oo 



lim e 



iujK 



U(-L, -oo)e 



-iujK 



-iujK 



= U (+oo,0) 
= Ub(0,-oo) 



(12) 



where Uq is the same as U, but with the self-interactions only (since these two factors 
correspond to the evolution of the projectile while outside of the target field). For 
the factor U(+L, —L), the change e"x _ — > x~ gives 



lim e 

id— >+O0 



iujK 



U(+L,-L)e- tL 



iuK 



exp 



ig / d 2 x ± x(x±)p(x ± ) 



(13) 



with 



dx + A (.t + ,0,£c_l) 



tt(oo) 



(14) 



p(x±) = J dx~ J + (0, x~ , x±) . 
Thus, the high-energy limit of the scattering amplitude is 

(An|^o(+co,0) exp ig J d 2 x± x(x±)p(x±) U (0, -oo)|a in ) . (15) 

A few remarks are in order at this point 

• Only the A~ component of the vector potential matters. 

• The self-interactions and the interactions with the external potential are 
factorized into three separate factors - this is a generic property of high 
energy scattering. 

• This is an exact result in the limit co — > +oo. 
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Eq. (15 1 is an operator formula that still contains the self-interactions of the 
fields to all orders. In order to evaluate it, one must insert the identity operator 
written as a sum over a complete set of states on each side of the exponential, 



4a ) =E(^l U 0(+ OO '°)l^) 
7,<5 



The factor 



(7in|exp ig j x(xj_)p(xj_) \S in )(5 in \U(0, -oo)|a in ) . (16) 
5^|*in>(*ia|tf(0,-oo)|a io ) (17) 



is the Fock expansion of the initial state: it accounts for the fact that the state a 
prepared at x + = — oo may have fluctuated into another state 5 before it interacts 
with the external potential. The matrix elements of Ug that appear in this expansion 
can be calculated perturbatively to any desired order. There is a similar factor for 
the final state evolution. 

The interactions with the target field are contained in the central factor, 
(7in| ex P •■■|^in)- In order to rewrite it into a more intuitive form, let us first rewrite 
the operator p as 

P a ^ =t *ij ^|^|^{^^ + ^-^-^ + '^^>^^ ) ^ 

-4 n (p + ,P±li)din(p + ,q±;j)e- i ^- q ^} , (18) 

where the tfj are the generators of the fundamental representation of the SU(N) 
algebra and the b, d, , d) the annihilation and creation operators for quarks and 
antiquarks. p a also receives a contribution from gluons, not written here, obtained 
with the generators in the adjoint representation and the annihilation and creation 
operators for gluons instead. This formula captures the essence of eikonal scattering: 

• Each annihilation operator has a matching creation operator - therefore, 
the number of partons in the state does not change during the scattering, 
nor their flavor. 

• The p + component of the momenta are not affected by the scattering. 

• Only the colors and transverse momenta of the partons can change during 
the scattering. 

Scattering amplitudes in the eikonal limit take a very simple form if one trades 
transverse momentum for a transverse position, by doing a Fourier transform. For 
each intermediate state (5 in | = fejx}| > define the corresponding light-cone 

wave function by : 

*5c({kt,Xi±}) = H[ d ^e- ik *^(5 in \U(0,-oc)\a in ) , (19) 
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where the index i runs over all the partons of the state <5. Then, each charged 
particle going through the external field acquires a phase proportional to its charge 
(antiparticles get an opposite phase), 

V S a({k?,Xix}) — > aJij.})JJZ7i(a!x) 

i 

Ui(x±) =Texp ig J dx + A~(x + , 0, Six)*"] . (20) 

(U is replaced by W for antiquaries, and t a by an adjoint generator for gluons.) The 
factors U in this formula are called Wilson lines. 



2.2. Dipole scattering at Leading and Next to Leading Order 

Let us now assume that the initial and final states a and /? both contain a virtual 
photon, as would be the case in DIS at leading order. The simplest Fock state that 
contributes to its wave function is a qq pair in a color singlet state, and the bare 
scattering amplitude can be written as 



« * g> * (x ± , y± ) * <°> (x ± , V± ) C/, fe (a x ) E/£ (y ± ) 

*(°)( :C x,yx) 2 tr[^x)^ t (y±)] ■ (21) 



The color trace arises because the wave function of the photon is diagonal in color, 
^ij ~ " ^ n diagram in the left hand side, the gray band represents the Lorentz 
contracted target field. 

It turns out that one-loop corrections due to the emission of a gluon inside the 
dipole are enhanced by logarithms of the dipole longitudinal momentum p + , leading 
to possibly large corrections proportional to a s log(p + ). In order to compute the 
dipole amplitude at next to leading order, we need to evaluate the following graphs 

<$> <$> 



> <T> <ZL> 



+ h.c. 



We will call real the terms where the gluon traverses the target field, and virtual 
those where the gluon is a correction inside the wavefunction of the incoming or out- 
going dipole. In the light-cone gauge A + = 0, the emission of a gluon of momentum 
k and polarization A by a quark can be written as 



% = 2 .< ^TT^ - (22) 

where it is assumed that the gluon is soft compared to the quark and where 
is the polarization vector of the gluon. Trading transverse momentum in favor of 
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transverse position, this vertex reads 



(27T)* 



'- 2 " 2tt (a^ - z x ) 2 



k\ 



(23) 



The other rule we need to compute these graphs is that when connecting the gluons 
to form the loop, one must sum over their polarizations, leading to a factor 



(24) 



Let us start with the virtual contributions. For instance, one gets 

\¥°\x x ,y x )\\r [tH a U{x x )tf{y x )] 

d 2 z x (xj_ - z±) ■ (x± - z x ) 
(2tt)2 (x x -z x ) 2 (x x -z x f 



x — 2a. 



dk+ 
~k+ 



(25) 



and 



<z[> 



& \x x ,y x ) tv[t a U(x x )tf(y x )t a ] 

d 2 z x (x± - z±) ■ (y ± - z x ) 



x4a 



dk+ 

IF 



(2tt)2 {x x -z x f{y x -z x f 



(26) 



Similar expressions can be obtained for the other virtual terms, and the sum of all 
virtual corrections is 

. _ ... 

d 2 z x 



C F a s 



dk^ 



(x x - y x ) 



(x x - z x ) 2 (y ± - z x f 



^°\x x ,y x ) ti[U{x x )U\y x j\ , 



(27) 



where we have introduced the Casimir of the fundamental representation of SU (N), 
C F = t a t a = (N 2 — 1)/2N. This (incomplete, since the real terms are still missing) 
result exhibits two pathologies : 

• The integral over k + is divergent. It should have an upper bound at p + , 
the longitudinal momentum of the quark or antiquark, 



= log(p+) = Y 



(28) 



When the rapidity Y is large, a s Y may not be small. By differentiating 
with respect to Y, we will obtain an evolution equation in Y whose solution 
resums all the powers (a s Y) n . 
• The integral over z x is divergent near z x = x x or y ± . This is a collinear 
singularity, due to emission of the gluon at zero angle with respect to the 
quark or antiquark. We will see shortly that it is cancelled when we combine 
the virtual and real contributions. 
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Let us now turn to the real terms. For instance, one has 
* (0) (*x,I/_l) 'tr [t a U(x ± )t b tf(y ± )} 
x4a 




~k+ 



— t Uab{z± \ x± -z ± nx ± -zjv ' (29) 



where U a b is a Wilson line in the adjoint representation, describing the eikonal 
scattering of the gluon off the target field. In order to simplify the color structure, 
first recall that 



t a U ab (z ± ) = U(z ± )t b UHz ± ) 
and then use the SU (N) Fierz identity, 



(30) 
(31) 



Thanks to these identities, the Wilson lines can be rearranged into 

tr [t a U(x ± )t b tf(y ± )] U ab ( Zl _) = *tr [U*(z±)U(x±)] tr [U(z±)U*(y ± )] 

-^*Wx ± )rf(v ± )] . (32) 

Interestingly, the subleading term in the number of colors (i.e. in 1/2N C ) cancels 
exactly against a similar term in the virtual contribution. Moreover, when we sum 
all the real terms, we generate the same kernel as in the virtual terms, 

(xj_ -y ± ) 2 
(xj_ - zj_) 2 (y ± - zj_) 2 ' 

In order to make the following equations more compact, let us introduce 

S(x ± ,y ± )^^-tr[U(x ± )U\y ± )] . 

In terms of this object, the sum of all the NLO contributions is 

(x± -y ± ) 2 



(33) 



(34) 



a s N 2 Y 
2ir 2 



* (0) (*l,V±)| 2 / d 



(xj_ - zj_) 2 (yj_ - z ± ) 2 
{S(xx,y ± ) - S(x±,z±)S(z±,y x )\ 



(35) 



while the bare scattering amplitude was |^°)(cc_i_, y±)\ N c S(x±,yj_) . By com- 
paring these two formulas, we conclude that 



dS(x ± ,y 1 _) 
dY 



2tt 2 
X ■ 



d 2 z ± 



(x± - y ± ) 



(x ± - z ± ) 2 (yj_ - z ± ) 2 
- S(x ± ,z±)S(z±,y ± )^ 
Note that, since S(x±, x±) = 1, the integral over z± is now regular. 



(36) 
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2.3. BFKL equation 



In the derivation of eq. ( 36 1 , we have not made any assumption about whether 
the amplitude S is large or small. Let us now focus on the dilute regime, which 
corresponds to small scattering amplitudes T(x±, y ± ) = 1 — S(x±, y±). By rewrit- 



ing eq. (36 1 in terms of T and by keeping only the linear terms, we obtain the 



Balitsky-Fadin-Kuraev-Lipatov equation (BFKL) PUGS 

dT(x±,y x ) = a s N c f 2 (x ± -y ± ) 2 

dY 2tt 2 J ± (xi - z x ) 2 (y ± - 2i) 2 



:{T(x x ,z x )+T(z ± ,y ± )-T(x ± ,y ± )} . (37) 



The BFKL equation has a fixed point, T = 0, but it turns out to be unstable. 
Moreover, one can see that it leads to violations of unitarity : 

• The mapping T — » a s N c J ■ ■ T has a positive eigenvalue w. 

• Therefore, generic solutions of the BFKL equation grow exponentially with- 
out bound as exp(ujY) when Y — > +oo. This leads to a violation of unitarity, 
since T-matrix elements should be bounded. 

This growth of the scattering amplitude with rapidity can be related to the increase 
at small x of the gluon distribution. Indeed, if one is still in the dilute regime, the 
forward scattering amplitude between a small dipole and a target made of gluons 
is proportional to 

T(x±,y ± ) cx \ x± - yj 2 xG(x, \x± - yj" 2 ) (38) 
where Y = ln(l/x). Therefore, the exponential behavior of T translates into 

x G{x,Q 2 ) ~ ^ ■ ( 39 ) 
x M 

2.4. Gluon saturation and BK equation 



Interestingly, the original equation (36), did not have this unitarity problem. When 
written in terms of T, it reads 

dT(x±,y ± ) _a s N c f ^ (x ± - y ± ) 2 



dY 2^ J - (x ± -z ± )*(y ± -z ± )2 

x {t(x_l, z±) + T{z ± ,y ± ) - T(x ± ,y ± ) - T{x ± , z±)T(z±, y ± )} , (40) 

an equation known as the Balitsky-Kovchegov equation (BK) 112IE2I One can see that 
it has two fixed points, T = and T = 1. The appearance of the second fixed point 
at T = 1 is due to the nonlinear term in T, that we had neglected in deriving the 
BFKL equation. Moreover, this new fixed point turns out to be stable. Therefore, 
the typical behavior of solutions of the equations it that T starts at a small value, 
grows exponentially with Y, and then saturates at T = 1. 
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The physical interpretation of what makes the solutions of the linearized BFKL 
equation grow exponentially, and why the solutions are bounded when one keeps 
the nonlinear term in T is rather transparent from the point of view of the gluon 
distribution in the target. This is illustrated in the left part of the figure IU At small 




Fig. 8. Left: gluon cascades from a valence quark in a hadron. Right: saturation domain as a 
function of the value of x and the atomic number A of the nucleusPB 



Y, i.e. in the dilute regime where T is small, one probes the valence quarks only. 
At larger values of Y, the dipole can interact not only with the valence quarks, 
but also with gluons radiated from the valence quarks. At even larger values of Y , 
these gluons themselves emit more gluons, in ladder-like diagrams. In this regime, 
the growth of the number of gluons is exponential in the rapidity Y. This very fast 
growth continues until the gluon density becomes too large. Roughly speaking, when 
the areal gluon density in the transverse plane multiplied by the cross-section for 
the recombination of two gluons by the process gg — > g becomes of order one, then 
the nonlinearities become important. Their effect is to tame the growth of the gluon 
density, so that the scattering amplitude T remains bounded. This phenomenon is 
known as gluon saturation. 

The criterion for the onset of gluon saturation can be expressed in terms of the 
resolution scale Q at which the gluons are observed. It reads 

OsQ^ x A- 2 / 3 xG(x,Q 2 ) > 1 , (41) 

tfgg-^g surface density 

where A is the atomic number in the case the target is a nucleus (otherwise it is 
unity). This inequality can be rearranged into 

Q 2 <Q 2 S = a ° xG ^f^ „ ^Vs^-o.s (42) 

saturation momentum 

The right hand side of this inequality is called the saturation momentum. The 
effects of gluon saturation are important below this scale, and are therefore more 
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visible when this scale is large. It is easy to estimate the A and x dependence of the 
saturation momentum. Because Q 2 is the ratio of a quantity proportional to the 
volume of the target by an area, it should scale like A 1 / 3 . Its x dependence can be 
inferred from the behavior of the gluon distribution in nucleons at small x: it grows 
like a power of 1/x, the numerical value ~ 0.3 of the exponent being extracted from 
HERA data for instance. The value of the saturation momentum as a function of 
x and A is shown in the right plot of the figure [8j and the dense/saturated regime 
comprises the volume under the surface. 

Note that the eq. (40) derived here resums only the leading log terms, since it 
is based only on the kernel at one loop. At this order, there is no running of the 
coupling constant. This effect arises via some parts of the two-loop kernel, that 
involve the beta function. These corrections are known by now^^U and play an 
important role in reaching a good agreement with experimental dataP^HSll 



2.5. Target average 

Until now, we have treated the target as a given patch of color field, produced by 
the constituents of the target, that are static over the timescale of the collision due 
to Lorentz time dilation. However, this picture implies that details of this color 
field depend of the precise configuration (position, color, ...) of the constituents of 
the target at the time of the collisions, that of course we do not know. Since this 
configuration changes from collision to collision, so does the target field. Therefore, 
instead of a single target field, one should instead have in mind an ensemble of 
such fields, corresponding to all the possible configurations of the constituents of 
the target. This means that all the quantities we have so far evaluated in a given 
target field should in fact be averaged over this ensemble of target fields, 

T -> (T) . (43) 

When performing this average at the level of the BK equation, we naturally get 

d (T{x ± ,y ± )) = a s N c f 2 (x± - Vl _) 2 

dY 2n 2 J z± {x±-z 1 _Y{y x -z x y 

x {(T(x ± , z ± )) + (T(z ± ,y ± )) - (T(x±,yj_)) - (T(x ± ,z ± )T(z ± ,y ± ))} .(44) 

The trouble with this equation is that it is no longer a closed equation, since the 
evolution with Y of the quantity (T) depends on the value of a new quantity (TT) . 
By the same method used to derive the evolution equation for (T) , we could derive 
an evolution equation for (TT), but its right hand side would contain the average 
value of an object containing three factors T, and so on. Therefore, instead of a single 
closed equation, we have now an infinite hierarchy of nested equations, known as 
the Balitsky hierarchy^ 

There is an approximation of the above hierarchy of equations, that leads us 
back to the original Balitsky-Kovchegov equation. This approximation amounts to 
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assuming that the average of a product of two dipole operators factorizes into the 
product of their averages, 

(TT) « (T) (T) . (45) 

Obviously, if this is true, then the first equation of the hierarchy is nothing but the 
BK equation. This mean-field approximation is believed to be valid for large targets 
such as nuclei, and in the limit of a large number of color. Although subject to this 
approximation, the BK equation is widely used in phenomenological applications 
because of its simplicity. 



2.6. Geometrical scaling 

The BK equation leads to an interesting property of the DIS total cross-section, 
called geometrical scaling. Firstly, let us write the 7*-target cross-section as follows, 



cr 7 » T = J dz J d 2 r J _\ip(q\z,r A _)\ 2 a d i P o\c(rj_} 



(46) 



with 



CTdipoicK) = 2 J d 2 X ± T(X X + y,Xi - —) . (47) 

ip(q\z, r±) is the light-cone wavefunction for a virtual photon of momentum q split- 
ting into a qq pair or transverse size r j_ , where the quarks carries the fraction z of 
the longitudinal momentum of the photon. This object is a purely electromagnetic 
quantity (at least at leading order), and is not interesting from the point of view 
of strong interactions, that are all contained in the dipole cross-section defined in 



cq. ( 47 1 . This equation defines the total cross-section between the qq pair and the 
target, expressed in terms of the forward scattering amplitude thanks to the optical 
theorem. Thus, the BK equation tells us about the rapidity dependence (i.e. the x 
dependence since Y = log(l/x)) of the DIS cross-section. 

For the sake of the argument, let us assume translation and rotation invariancc, 
and define 

N(k ± ) = 2tt / d 2 Xl _ e ik ^ r 7 ^ ^)) . ( 48 ) 



From the Balitsky-Kovchegov equation for (T), we obtain the following equation 
for N 

dNjkQ = a s N c 
dY TT 



x(-9 L )iV(fc ± )-7V 2 (fc ± ) , (49) 



where we denote 

L = Mk 2 /kl) 



X (t) = 2V»(1) - ^(7) - V>(1 - 7) 
d\nT(z) 

ip(z) = — . (50) 

dz 
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(T(z) is Euler's Gamma function.) The form (49) of the BK equation is particularly 



simple, because its nonlinear term is just the square of the function N. 



The function x{l) nas a minimum at 7 = 1/2. Eq. (49) can be approximated by 
expanding x{l) to quadratic order around the minimum. This amounts to keeping 
only derivatives with respect to L up to second order, i.e. to a diffusion approxi- 
mation. The physical content of this equation becomes more transparent if one also 
introduces the following variables, 



t ~ Y 



L+^ X "(l/2)Y, (51) 



in terms of which eq. ( 49 ) becomes^^H 

d t N = d 2 N + N-N 2 . (52) 

This equation is known as the Fisher-Kolmogorov-Petrov-Piscounov (FKPP) equa- 
tion, and it describes the physics of reaction-diffusion processes. These are processes 
involving a certain entity A that can do the following: 

• An entity A can hop from a location to neighboring locations. This is de- 



scribed by the diffusion term in the right hand side of eq. ( 52 ) . 

• An entity A can split into two, increasing the population by one unit. This 
is described by the +N term. 

• Two of these entities can merge into a single one, reducing their population 
by one unit. This is described by the term —N 2 . 

The FKPP equation has two fixed points, N = which is unstable, and N = 1 
which is a stable fixed point. Note that the loss term, —N 2 , is essential for the 
existence of this stable fixed point. Moreover, for rather generic initial conditions, 
the solutions of the FKPP equation are known to behave like traveling waves at 
asymptotic times 

N(t,z) ~ N(z-2t). (53) 

t— >+oo 

In other words, these solutions propagate in the +z direction at the constant velocity 
dz/dt = 2, and instead of depending separately on z and t, they depend only on the 
combination z — 2t. Going back to the dipole scattering amplitude T, we see that it 
does not depend separately on Y and on the dipole size r± , but on the combination 
Q s (Y)r±, where Q S (Y) has the following Y dependence 

Ql(Y) = k 2 y-5irN) e ^x"im-y)Y . (54) 

This quantity is nothing but the saturation momentum introduced by a simple 
semi-quantitative argument earlier. Here, we see it emerge dynamically from the 
evolution equation itself. This is a major feature of gluon saturation: the ability 
to generate a dimensionful scale is a consequence of the nonlinear aspect of the 
saturation phenomenon. Moreover, thanks to the fact that this scale increases with 
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Y, i.e. with the energy scale, one may hope that gluon saturation falls in the realm 
of weakly coupled physics at sufficiently high energy. 

And at the level of the 7*-target cross-section, this implies that the cross-section 
depends only on Q 2 /Q 2 (Y) : rather than Q 2 and Y separately (this is strictly true 
assuming that we neglect any other dimensionful parameters in the problem, such as 
the quark masses that enter in the virtual photon wavefunction ip - this is legitimate 
for the light u, d and s quarks). This scaling property of the DIS total cross-section 
has been observed in experimental data, as shown in the right plot of the figure 
|9j What makes gluon saturation a satisfying explanation of this phenomenon is 




Fig. 9. DIS data at HERA. Left: displayed conventionally as a function of x and Q 2 . Right: 
displayed as a function of r = Q 2 /Q 2 a (x)EMB 



that geometrical scaling is a very robust consequence of evolution equations with a 
nonlinear term that tames the growth of the amplitude, since any realistic initial 
condition falls leads to traveling wave solutions. Interestingly, the scaling shown in 
the right plot of the figure [9] also works for DIS off nucleip^ provided one scales the 
saturation momentum by the appropriate factor A 1 / 3 (see the NMC data points). 
For large nuclei like the ones employed in heavy ion collisions (Gold at RHIC and 
Lead at the LHC), the nuclear enhancement of the saturation momentum is A 1 / 3 ps 
6. From fits to existing data on deep inelastic scattering on various targets, one can 
get a more precise idea of the numerical value of the saturation momentumEUH 
(this is how Q s has been obtained in the figure [8]). 
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3. Introduction to the Color Glass Condensate 
3.1. Other elementary reactions 

In the previous section, we have studied Deep Inelastic Scattering and its energy 
evolution from the point of view of the projectile (in this case, the virtual photon). 
When going to NLO, the correction was applied to the qq dipole, while the target 
field was held fixed. In other words, all the energy evolution was applied to the 
dipole. This is perfectly legitimate: since cross-sections are Lorentz invariant objects, 
they do not depend on the frame. The kind of calculation we have performed in the 
previous section amounts to observing the reaction in a frame where the target is 
nearly at rest, and the boost to go to higher energy is applied entirely to the virtual 
photon. In the case of DIS, the interest in doing so is obvious: the target is very 
complicated (a nucleon or a nucleus), while the projectile is a rather elementary 
object (a virtual photon). Thus it seemed natural to look at the energy evolution 
by boosting the object we understand best. 

There is another situation where one can proceed in this way, that involves 
colliding a dilute projectile on a dense target. This is the case for instance in proton- 
nucleus collisions, at forward rapidities in the direction of the proton beam. In this 
kinematical configuration, one probes the dilute regime of the proton, and the dense 
regime of the nucleus. In this case, the proton can be describe as a dilute beam of 




Fig. 10. Examples of elementary processes involving a projectile quark or gluon and the field of 
a target nucleus. Each black dot represents the Wilson line associated with the eikonal scattering 
of the corresponding parton. 

quarks or gluons whose flux is given by the usual parton distributions, and the 
elementary reactions we must consider are the interactions between a quark or a 
gluon and the field of the target nucleus. Examples of processedHHUHUlflll that have 
been evaluated in this way are represented in the figure [T0| 

However, this approach suffers from a serious limitation, since it is practical 
only when the projectile is much simpler than the target. For symmetric collisions, 
such as nucleus-nucleus collisions. In this case, there is no advantage in treating 
differently the two nuclei. In fact, doing so is arguably rather unnatural. As we 
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shall see shortly, the Color Glass Condensate description of gluon saturation allows 
a perfectly symmetric description of the two nuclei in these collisions. 

3.2. Color sources and quantum fields 

Before we turn to the study of nucleus-nucleus collisions, let us consider again 
DIS, but this time from the point of view of the target, i.e. when we go to higher 
energy, we apply the boost to the target while leaving the projectile wavefunction 
unchanged. In the previous, we simply assumed that the target can be represented 
by a certain ensemble of color field configurations, that we did not specify. If we 
adopt the point of view that consists in applying the boost to the target, then this 
ensemble of fields must change with energy for the cross-section to have the correct 
energy dependence. 

At the basis of the Color Glass Condensate framework is the fact that the fast 
partons of a hadron or nucleus suffer Lorentz time dilation and therefore do not 
evolve during the short duration of a collision. They can thus be considered as time 
independent objects moving at the speed of light along the light-cone, and the only 
relevant information we need to know about them is the color charge they carry. 
This can be encoded in a color current of the form 

JZ = 5» + p a (x-,x ± ) . (55) 

This formula keeps only the longitudinal component of the current (the + compo- 
nent in light-cone coordinates), and neglects all the other components, because only 
this component is enhanced by the Lorentz factor. The lack of x + dependence stems 
from the assumption that these color charges are time independent. The function 
p a describes the spatial distribution of these color charges, both in x~ and in x±. 
Note that the support of its x~ dependence is very narrow, due to Lorentz contrac- 
tion. At leading order, the target field can be obtained by solving the Yang-Mills 
equation with the source J^, 

[D„,F^] = J". (56) 

In the Lorenz gauge, defined by d^A^ — 0, this equation can be solved analytically 
for the above current, and one obtains 

A~=A i =0 , A+(x) = \-p{x~,xx). (57) 

Therefore, specifying the ensemble of target fields is equivalent to specifying the 
ensemble of the functions p. Thus, the CGC framework must be supplemented by 
a probability distribution 

The ideas behind the Color Glass Condensate originate in the McLerran- 
Venugopalan modeP^^ that, in addition to the description of a dense object in 
terms of classical color sources, also argued that their distribution W[p] should be 
nearly Gaussian, with only local correlations, 

(p a (x-, Xl _)p b (y-,y x )) =p 2 5 ab 5{x- - y~)S{x ± - y ± ) . (58) 
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The justification of this model is that, in a highly Lorentz contracted nucleus, there is 
a large density of color charges at each impact parameter. The charges from different 
nucleons are uncorrelated, and thus by the central limit theorem, the resulting 
distribution should be approximately Gaussian for a large nucleus. The locality of 
the correlations is also a consequence of the fact that charges in different nucleus 
are not correlated. 




k 



Fig. 11. Left: illustration of the McLerran-Venugopalan model for a large nucleus. Right: sepa- 
ration in longitudinal momentum between fields and sources (here for a target moving in the — z 
direction - hence the cutoff on the fc~ component of the momentum). 



3.3. Cutoff dependence and renormalization group evolution 

The justification for treating the constituents of the target as static color sources is 
Lorentz time dilation. Arguably, this is only justified if these partons have a large 
enough longitudinal momentum. The partons that are too slow must be treated in 
terms of the usual gauge fields (if they are gluons). Therefore, the CGC should be 
viewed as an effective theory with a cutoff A + , such that 



Partons with fc + > A + are treated as static sources, via eq. (55 1 



• Partons with fc + < A + are treated as standard quantum fields. 

Because these two types of degrees of freedom have well separated longitudinal mo- 
menta, the interactions between them can be approximated by an eikonal coupling 
J^A 11 , and the action that describes the complete effective theory is 

S = -\f F tlv F^+ j J»A^ . (59) 



Yang-Mills fast partons 

In the CGC framework, the expectation value of an observable is obtained by first 
computing the observable for a fixed configuration p of the color sources of the 
target, and then by averaging it with the distribution 

(O) = j [Dp] W[p] 0[p] . (60) 

The cutoff A + that has been introduced to separate the fast and the slow partons 
is arbitrary, and observable quantities should not depend on it. At leading order, 
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this is always the case. But when one computes loop corrections, the cutoff enters 
in the calculation as an upper limit on the longitudinal momentum that runs in 
the loop, in order to prevent the quantum modes that are integrated in the loop 
to duplicate modes that are already introduced via the static source p a . In general, 
this leads to a logarithmic dependence of the loop corrections on the cutoff A+P^l 
The only way to have A-independent observables at the end of the day is to let 
the distribution W[p] depend itself on A, precisely in such a way that the two A 
dependences cancel, 

W[p]^W A [p\. (61) 

Moreover, it turns out that this can be achieved with the same distribution W A [p] 
for a wide range of different observables, with a A dependence governed by the 
so-called JIMWLK equationj£2H2£l of the form 



dW 



~nW A , (62) 



d log A 

where H (the JIMWLK Hamiltonian) is a quadratic operator in functional deriva- 
tives with respect to p. Note that the A dependence is often traded for a rapidity 



dependence, by using Y = log(A). The evolution equation (62) can be seen as a 
renormalization group equation, that describes how the color charge content of the 
target changes as one changes the longitudinal momentum A down to which they 
are considered. When the cutoff A is large, and close to the longitudinal momentum 
of the target, only its valence partons are included in the description in terms of 
p. Lowering A means that one includes in this effective description more and more 
soft modes (sea partons). 

The JIMWLK equation does not predict by itself what the distribution W A [p] 
is for a given target, only how it changes when on lowers the cutoff. It must be 
supplemented by an initial condition at some A in order to lead to definite results. 
This initial condition is by essence non-perturbative, and must be modelled^] When 
the evolution of the distribution W A [p] is considered, the status of the McLerran- 
Venugopalan model changes a bit, because the Gaussian distribution proposed by 
the MV model is not a fixed point of the JIMWLK equation. Instead, the MV model 
is often viewed as a plausible initial condition for a large nucleus at large cutoff (i.e. 
close to the valence region of the target). 

Although it is necessary to specify an initial condition in order to solve the 
JIMWLK equation, it should be stressed that its solutions have a universal scal- 
ing property when evolved far enough from the initial cutoff scale: the equation 
generates an intrinsic momentum scale, the saturation momentum Q s , and in this 
scaling regime all the correlation functions of Wilson lines depend on the transverse 
coordinates only via combinations such as Q s x±. This scaling of the solutions of the 

In this sense, the JIMWLK equation has exactly the same status as the DGLAP equation for the 
Q 2 dependence of the ordinary parton distributions. 
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JIMWLK equation is what led to asymptotic traveling wave solutions for the BK 
equation. This behavior of the solutions mean that they tend to forget the details 
of their initial condition, if evolved sufficiently far from the starting scale. This sug- 
gests that predictions of the CGC should be less sensitive to the model employed 
for the initial condition when applied to the study of collisions at very high energy. 



3.4. Relationship between the Balitsky's and CGC formulations 

The correspondence between Balitsky's hierarchy and the CGC formulation can be 
summarized by the following identity, 

(0) Y = J [Dp] W [p] Y [p] = J [Dp] W Y [p] O [p] • (63) 

v ' > v ' 

Balitsky's description CGC description 

In the Balitsky's approach, the rapidity dependence is obtained by applying a boost 
to the observable, while the distribution of sources is kept unchanged. In the CGC 
approach, the observable is a fixed functional of the sources, and the boost leads to 
a change in the distribution of the sources. In the simple case of collisions between 
an elementary projectile and a dense target, where observables can be expressed in 
terms of Wilson lines of the target field, the equivalence between these two points 
of view follows from the existence of a universal operator "H such that 

^^uLl] OM- (64) 



dY 



dp 



Moreover, this operator is self-adjoint, which means that one can "integrate by 
parts" and transpose its action from the observable to the distribution leading 
to the JIMWLK evolution equation. 



3.5. Example : DIS in the CGC framework 

As an illustration of the use of the CGC, let us reconsider Deep Inelastic Scat- 



tering in this framework (see the figure 12). We present here only a sketch of the 
calculation, but not the detailsEHEH At LO, the photon must fluctuate into a qq 
dipole in order to interact with the field of the target. As with the BK approach, 
the forward scattering amplitude of this dipole off the target can be expressed in 
terms of Wilson lines, 

T LO (x ±: y ± ) =l--L t r(U(x ± )U\y ± )) 

U(x ± ) = P expig J dz + A~(z + ,x±) , (65) 

and the main difference resides in the fact that the target field is not arbitrary but 
obtained by solving the Yang-Mills equation 

{V^^}=5»- p(x+,x ± ) . (66) 
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X 

Fig. 12. Deep inelastic scattering in the CGC framework. The green arrow represents the tra- 
jectory of the target, and the support of its color sources. Left: Leading Order: the virtual photon 
fluctuates into a qq dipole that scatters off the color field of the target. Right: Next to Leading 
Order: the photon fluctuates into a qqg state. 

Then the cross-section is obtained by averaging over all the configurations of the 
color source p, with the distribution At LO, the difference with the BK point 

of view is purely semantic, since we are simply trading the (arbitrary) distribution 
of the target fields for the (arbitrary) distribution of the color sources. 




fields »« — sources — » 

1 1 h 

a; a _p 

5T T 

NLO LO 



Fig. 13. Infinitesimal step in the cutoff evolution, from Aq to Ai. 



At NLO, one must now consider corrections such as the one represented in the 
right part of the figure 12 (there are several graphs, only one of them is represented 
there) . The longitudinal momentum of the gluon in the calculation of these correc- 
tions should be limited at the upper value A^~, to avoid over-counting modes that 
are already included via the sources p. In practice, it is more convenient to compute 
only the correction <5T NLO to the scattering amplitude due to the field modes in a 
small strip A^ < k~ < Aq\ The important part in this discussion are the terms 
that have logarithm divergences when A^ <C . One can show that they take the 
following form, 

5T Nh0 (x x ,y ± )=ln(^j U T LO (x ± ,y ± ) , (67) 

where H is an operator that acts on the p's. Besides the fact that they depend 
on an unphysical cutoff, these logarithms question the validity of the perturbative 
expansion. Indeed, the NLO corrections are suppressed by one power of the coupling 
constant a s compared to the LO result, but the appearance of these possibly large 
logarithms can compensate the smallness of the coupling. This is why one should 
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consider only a small slice of longitudinal momenta at a time. It turns out that 
these logarithms can be hidden by averaging over the color sources, thanks to the 
identity 



T h0 +6T^ ) A _={T h0 ) A _, (68) 
where we use the shorthand 

("■>A = J[Dp}W A [ P } ■•■ , (69) 
provided the distributions of sources at the scales Aq and Aj~ are related by 

W A - . (70) 



W A - 



I + lnf^M H 



Ar 



The meaning of eq. ( 68 1 is that the sum of the LO+NLO contributions in the 
original effective theory (that has its cutoff at the scale Aq ) is equivalent to the 
LO only in a new effective theory that has its cutoff at the lower scale A]", and 



a new distribution of sources given by eq. ( 70 1 . Note that eq. ( 70 1 is nothing but 
the JIMWLK equation for an infinitesimal change of the cutoff scale. Moreover, 
the right hand side of eq. ( 68 1 has no dependence on the former scale Aq~ . This 



means that in the left hand side, this dependence must cancel between the NLO 
contribution and the scale dependence of the distribution W[p\. 

The previous process, by which we integrated out the logarithms coming from 
a small slice of field modes, can be repeated indefinitely by considering a sequence 
of lower and lower cutoffs, 

■ • • < A" < • • • < Ar < Aq . (71) 

At each step, new logarithms are generated, that can be absorbed by defining a 
new effective theory at the lower scale A~ +1 and a new distribution of sources at 
this scale. The dependence on the cutoff at the previous scale A~ disappears in 
this process. We must repeat this procedure until we reach a cutoff scale which is 
lower than all the physically relevant longitudinal momentum scales: at this point, 
lowering the cutoff further only introduces new sources that are too slow to be 
relevant in the observable under consideration - the observable does not depend on 
these sources, and thus becomes independent of the cutoff] 

4. Factorization in high-energy hadronic collisions 
4.1. Introduction 

In the previous section, we have introduced the Color Glass Condensate effective 
theory, and illustrated it by considering again the example of Deep Inelastic Scat- 
tering. In this context, the CGC is a mere rephrasing of the physics that was already 



g In the case of DIS, the dependence on the cutoff disappears once it becomes smaller than xP , 
where P~ is the longitudinal momentum of the target and x = Q 2 /(2P ■ Q). 
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present in Balitsky's hierarchy. Roughly speaking, the CGC puts the emphasis on 
the energy evolution of the distribution that describes the target, rather than on 
radiative corrections to the operator that is being measured, but the end result is 
exactly the same. 

The situation is qualitatively different in the case of nucleus- nucleus collisions. 
In this case, it would be highly desirable to have a framework in which the projectile 
and target can be treated on the same footing, and the CGC appears as a good can- 
didate for achieving this. In the CGC framework, we expect that each of the nuclei 
will be described by static color sources p\ and pi-, moving in opposite directions 
along the light-cone, with their respective distributions and H^[p2]- 

In several worksp^H the CGC has been used as follows in order to compute 
observables in heavy ion collisions : 

• Pick randomly two color sources pi and P2 (one for each nucleus), for in- 
stance by using the Gaussian distributions encountered in the MV model. 

• Solve (numerically) the classical Yang-Mills equation with two sources, 



with an initial condition such that the color field vanishes at x 



(72) 



o 



-CO. 



• Evaluate the observable of interest on the classical color field obtained 
by solving the previous equation. For instance, the gluon spectrum is ob- 
tained as 



dYd 2 p ± 



• Repeat these steps in order to perform a Monte-Carlo average over the 
distributions for pi and p2- 

As an illustration, we show in the figure [14] some numerical results for the single 
inclusive spectrum of the gluons produced in a nucleus-nucleus collision, obtained 
by following this procedure. 




Fig. 14. Transverse momentum dependence of the single inclusive gluon spectrum.^] 
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The goal of this section is to provide a theoretical justification for the above 
procedure, and in particular to answer the following questions : 

• Does it correspond to the leading term in an expansion in powers of the 
coupling constant? 

• What is the justification for the use of retarded boundary conditions? 

• In the case of a collision between two nuclei, can one still absorb the loga- 
rithms that arise in loop corrections by letting the distributions and 
Wi\p2\ evolve according to the JIMWLK equation? 



4.2. Power counting 




Before we can start answering these questions, we need a way to assess the order 
of magnitude of a given graphlHEn] Typical graphs in a nucleus- nucleus collisions are 
made of several disconnected components (see the graph on the right of the figure 
[6]). In order to perform this power counting, it is sufficient to consider only one of 
these connected components, as illustrated in the figure [15} In this representation, 
the circular dots represent insertions of the sources p\ or p2 - In the saturated regime, 
their contribution to the power counting is g , since the gluon occupation number, 
proportional to p 2 , should be of order g~ 2 . In addition, three-gluon vertices bring 
one power of g and four-gluon vertices a factor g 2 . 

With these rules, one obtains the following order of magnitude for a connected 



graph such as the one represented in the figure 15 

\g n - +2n - , (73) 
where n E is the number of external gluons (the gluon lines terminated by an arrow 



in the figure 15 ) and n L the number of independent loops in the graph. 

Interestingly, the order of the graph does not depend on the number of internal 
lines and on the number of sources pi.2 attached to the graph. The latter property 
is specific of the dense regime, where the sources are of order g^ 1 . Indeed, each 
insertion of a source necessarily brings a factor g for the vertex at which this source 
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is attached to the graph; this vertex cancels the g _1 coming from the source and 
therefore inserting a source in the graph does not cost any power of g. This means 
that, for graphs with a fixed number of external lines and a fixed number of loops, 
there is an infinity of topologies of the same order in g that differ in the number of 
sources they contain. Thus, computing any observable in the CGC framework, even 
at Leading Order, requires that one resums this infinite series of terms. 
For instance, the inclusive gluon spectrum has the following expansion 

dNt 1 



d 3 p g 2 



co + ci g 2 + c 2 g 4 - 



(74) 



where the coefficients cq,c\,--- are themselves series that resum all orders in 
(gp 13 ) n . For instance, 

oo 

CO = E c °> n (S , / J i,= ) n ■ ( 75 ) 
n=0 

Our goal in the rest of this section is to calculate the complete zeroth order term, 
co/g 2 , and a subset of the higher order terms. 



4.3. Bookkeeping 



Among the observables that one could possibly consider, the inclusive observables 
are especially simple because they do not veto any final state. This is the case for 
instance of the single inclusive gluon spectrumPJ defined in terms of the transition 
amplitudes as 



d 3 p 



oo 

E 

n=0 



(n + 1) 



(n + l)\ 



n part, phase-space 



(pPl ■ "PnoutlOin) 



(76) 



where we use the shorthand c£<£> = d 3 p/(2ir) 3 2p to denote the invariant phase-space 
of a final state particle. 

In order to sum all the relevant graphs at a given order in g 2 , we need tools 
to list and manipulate them. For the sake of th discussion in this section, we will 
consider observables related to particle spectra, such as the inclusive gluon spec- 
trum, correlations, etc... All these observables can be obtained from a generating 
functional 

1 



m - E ■ 



d$i • • • d$ r 



z(Pl) ■ ■ ■ z(P n ) <Pl"-P„out|0in) 



(77) 



h By using the completeness of the final states, one can check that this spectrum is also the 
expectation value of the number operator, 

dN\ . I | I i 

(Oin|ai u t(p)<W(p)|O in ) . 
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obtained by weighting each final state particle by an arbitrary function z(p). The 
diagrammatic interpretation of this generating functional is illustrated in the fig- 
ure 16 The product of an amplitude by its complex conjugate can be interpreted 





Fig. 16. Left: product of a term of the amplitude with a term from the complex conjugate 
amplitude. Right: interpretation of this product as a cut vacuum graph. 



as performing a cut through a vacuum graph (i.e. a graph without any external 
gluon) EHHIl To construct the generating functional F[z], one simply needs to weight 
each cut propagator (the propagators with a cross in the right part of the figure 
16) by the appropriate function z(p). From this generating functional, the single 
inclusive gluon spectrum introduced above is given by 

d/Vi _ SF[z] 
d 3 p 



8 z(p) 



(78) 



z=l 



In fact, quite generally, inclusive observables are derivatives of F[z] at z = 1, while 
exclusive observables are obtained via derivatives at other values of z. 

A crucial property in the CGC effective theory is unitarity, i.e. at a basic level 
the fact that the sum of the probabilities for all the possible final states in a collision 
should equal unity. In terms of the generating functional F[z], unitarity is encoded 
in a very simple way, F[l] = 1. The fact that inclusive observables are given by 
derivatives at z = 1, combined to this property due to unitarity, is the reason why 
they are much simpler than exclusive observables. In more physical terms, inclusive 
observables are simpler because many graphs cancel thanks to unitarity. 

By using reduction formulas for the transition amplitudes and their complex 
conjugate, one can write the generating functional F[z] as follows 

F[z]=e c M Z\j + }Z*\j_}\ j+=j _ =j , (79) 

where Z[j] is the usual generating functional for the time-ordered Green's functions 
that enter in the reduction formulas,^ and where 

__S S_ 

Jy Sj+ix) 5j-iy) 

G° + _ix,y) = J ^ e*-(-») zip) 2 7 r0(-p°)% 2 ) ■ 

(Note that here we have simplified the notations by not writing color, Lorentz 
and polarization indices - when dealing with gluons, these indices are of course 



CM = 



Jd 4 xd 4 y G° + _(x,y) U X U % 
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necessary.) Eq. (79), albeit being rather formal, is useful to obtain more explicit 



expressions for observables than can be derived from the generating functional. At 
this point, it is crucial to note that 

e c W Z\j + ]Z*[j.]=Z\j + J_], (80) 

where Z\j+,j-] is the generating functional for path-ordered Green's functions in 
the Schwinger-Keldysh formalism. 94 95 





Fig. 17. Diagrammatic representation of the two terms in the right hand side of eq. \81\ . Left: 
term in A+A_. Right: term in The crossed wavy line represents the on-shell momentum p. 



For instance, one can use eq. (79 1 in order to obtain the following expression for 
the single inclusive spectrum, 

dNi _ SF[z] 



= / d A xd A ye^ x -^ U x U y A + (x)A_(y) + G + -(x,y) , (81) 



d 3 p S(p) 

where A± and G_j are respectively the one-point and two-point Green's functions 

in the Schwinger-Keldysh formalism. The two terms in this formula are illustrated 
in the figure [17] This formula is true to all orders, and evaluating the spectrum at a 
given order amounts to evaluated these Green's functions to the desired order (the 



power counting rule of eq. ( 73 ) can be applied to these functions, respectively with 
n E = 1 and n E = 2). 



4.4. Leading order and retarded classical fields 



When we apply the power counting rule ( 73 1 to A± and G^ , we see that they have 

the following expansion in g 1 



A+ = 



G, 



a 
9 

h - 



+ aig + a 2 g + 



h 9 



">2 9 



(82) 



where the and hi are coefficients of order one. Therefore, to evaluate the spectrum 
at leading order, we need only the LO of A±, and we can disregard completely 



the term with G_| , as illustrated in the left figure 18 Despite this simplification, 

one still needs to sum an infinite set of graphs. Indeed, all the tree topologies 
contribute to A± at the order 1/g, and for all the internal vertices, we must sum 
over all the possible assignments of the types + or — , per the rules of the Schwinger- 
Keldysh formalism. 94 95 Since the color sources are identical on the two branches 
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Fig. 18. Left: only contribution to the single inclusive spectrum at leading order. Right: Aj r {x) 
at leading order. 



of the Schwinger-Keldysh closed time contour, this sum only gives the following 
combinations of propagators 

G\ + { P ) G° + _(p) = G\(p) G a __(p)= G°» , (83) 

i.e. the retarded propagator. Thus, we only need to sum all the tree diagrams built 
with retarded propagators, as illustrated in the figure [19] The result of this sum of 

- + K ^ 4^ 

Fig. 19. Beginning of the series of tree diagrams that contribute to A± at LO (only cubic vertices 
are included here, but of course in QCD one would also have four-gluon vertices). The black dots 
denote the color sources pi or p2, and the arrows on the lines indicate that the propagators are 
retarded. 



tree graphs is well known: it is the solution of the classical field equations of motion 
that vanishes in the remote past. In the case of QCD, these are the Yang-Mills 
equations, 

[V ll ,T» v ]=5 v -px+5 v +p 2 , lim A»{x) = . (84) 

This result is what justifies the procedure described in the introduction of this 
section in order to compute the inclusive gluon spectrum. 

Interestingly, at LO, the n-gluon inclusive spectrum factorizes into products of 
single gluon spectra, 



d 3 p 1 ■ ■ ■ d 3 p n 



dm 



d 3 p 1 



dN x 



d 3 P r , 



(85) 



as can be checked by simple power counting arguments. 



4.5. Next to Leading Order, logarithms and factorization 

The previous discussion justifies partly the procedure employed to compute observ- 
ables in the CGC framework. But since it is based solely on LO considerations, it 
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says nothing about the logarithmic cutoff dependence that may arise in loop cor- 



rections. Indeed, in eq. (74), it is implicitly assumed that the coefficients cq,ci, ■ ■ ■ 
are of order one. But these coefficients in fact contain logarithms of the cutoff, and 
should be more accurately written as 

ci = Cio + en InA* 

c 2 = c 20 + c 2 i log A* + c 2 2 log 2 A* 

Leading Log terms 

At each order, there can be powers of logarithms up to the number of loops. The 
terms that maximize the number of logarithms are called the Leading Log terms, 
and are the most important. 

On a more phenomenological side, the above LO result for the gluon spectrum 
is insufficient because at this order the spectrum has no dependence on the rapidity 
of the produced gluon. The two issues are in fact closely related : logarithms of the 
cutoff appear at NLO, and in order to cancel them one must let the distributions 
of sources Wif/Oi] and Wi\p-i\ evolve according to the JIMWLK equation. It is this 
evolution that gives the spectrum its rapidity dependence. 




Fig. 20. Causality argument explaining why the distributions W[p] are universal in nucleus- 
nucleus collisions. 



Let us first start with a qualitative argument to explain why it should be possible 
to absorb the logarithms of the cutoff into universal distributions W[p] in the case 
of nucleus-nucleus collisions (in fact the same distributions as those encountered in 



DIS). This argument, based on causality, is illustrated in the figure 20 and goes as 
follows 

• The duration of the collision is very short, and decreases as the inverse of 
the collision energy. 

• In contrast, the radiation of the soft gluons responsible for the logarithms 
takes a much longer time. Therefore, this radiation cannot happen during 
the collision itself - these gluons must be emitted before the collision. 
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• Before the collision, the two nuclei are not in causal contact. Therefore, 
what happens inside the first nucleus cannot be influenced by the second 
nucleus, and conversely, ft should also be independent of the quantities that 
an observer may measure in the final state (provided that this measurement 
does not lead to discarding some class of events - hence the special role of 
inclusive observables). Therefore, each nucleus should be described by a 
universal distribution W[p]. 




Fig. 21. Left: diagrams in the gluon spectrum at LO. Middle and Right: diagrams in the gluon 
spectrum at NLO. Here, the graphs are represented in coordinate space (with the two light-cones 
figuring the trajectories of the colliding nuclei). The green and red dots are the color sources carried 
by the two nuclei. 



Now, let us consider the NLO corrections to the single inclusive gluon spectrum, 
in order to see more precisely the structure of the logarithms. There are two types 
of contributions at NLO, represented in the middle and right graphs of the figure 



21 (for comparison, we have represented in the left graph a typical contribution at 

LO). The graph in the middle involves the two-point Green's function (x,y), 

that makes a first appearance at NLO, while the graph on the right contains a 
one-loop correction to the factor A+(x) (there is a third class of graphs, not shown 



in the figure 21 that contain a one-loop correction to A-(y)). 

The calculation of the NLO correction to the gluon spectrum is rather in- 
volved and will not be detailed here. Instead, we just sketch the main steps in 
this study™] 

i. Express the single gluon spectrum at LO and NLO in terms of classical 
fields and small field fluctuations. A crucial property is that these objects 
all obey retarded boundary conditions, 
ii. Write the NLO terms as a perturbation of the initial value of the classical 
fields. If Af i (x) is a solution of the classical Yang-Mills equations, and a^(x) 
a small (in the sense that <§C A^, so that the equation of motion of a M can 
be linearized) perturbation around this solution, they are formally related 
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by 



a" (a;) 



J [a(u)T u ] A»(z), 



(86) 



where T u is the generator of the shifts of the initial value of A 1 on some 
surface E (roughly speaking, one can see T u as a derivative with respect 
to the initial value of the field at the point u G E). Thanks to this identity, 
one can obtain the following formula for the gluon spectrum at NLO 









d 3 p 


NLO 





G(u, v)T u T v 



/3(u)T t 



dm 
d 3 p 



(87) 



where f3 and £7 are one-point and two-point functions that can be com- 
puted analytically. The main interest of this formula is that it provides a 
factorization in time of the NLO spectrum: the operator inside the square 
brackets depends only on the fields under the surface E, while the fac- 
tor dNi/d 3 p\- LO on which it acts depends only on what happens above E. 
Therefore, by choosing appropriately the surface E, this formula separates 
what happens in the two nuclei before the collision from the collision itself. 
Note also that in this formula, the first factor can be calculated analyti- 
cally, while the second cannot because it involves solutions of the classical 
Yang-Mills equations that are not known analytically. 




Fig. 22. Choice of the surface £ to extract the initial state logarithms. 



Choose E as in the figure 22 When u, v are both on the t = z branch of 
the light-cone, one has : 

i J g(u,v)T u T v + J (3(u)T u = logA+ x Hi + finite terms , 



where T-L\ is the JIMWLK Hamiltonian of the nucleus that moves in the +z 
direction. If the two points u,v are on the other branch of E, we obtain a 
logarithm of A~ , whose coefficient is the JIMWLK Hamiltonian H 2 of the 
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other nucleus. If the points u, v are on different branches of the light-cone, 
there is no logarithnj^J Therefore, at leading log accuracy, we have 



IV. 



dNi 
d 3 p 



Leading Log 



log (A+) Hi + log (A-)ft s 



d 3 p 



(89) 



Like in the case of DIS, the logarithms can be hidden by integrating over 
the sources p\ and p2, with distributions Wi[/9i] and Wi\p2\ that obey the 
JIMWLK equation, with the JIMWLK Hamiltonian of the corresponding 
nucleus. To see this, one simply uses the self-adjointness of the JIMWLK 
Hamiltonian in order to transpose its action from the observable on the 
distributions W [/£>]. The final formula at Leading Log for the gluon spectrum 
is 



d 3 p Leading Log 



[Dp, Dp 2 ] Wi[ Pl ] W 2 [p, 




(90) 



The factorization formula (90 1 ends the justification of the procedure employed in 



the literature in order to evaluate the gluon spectrum in the CGC framework. Sim- 
ilar factorized formulas also exist for n-gluon inclusive spectra, or for observables 
like the energy-momentum tensor. In each of these extensions, it is the same dis- 
tributions W[p] that enter in the formula, therefore providing examples of their 
universality. In order to extend this factorization to other observables, the central 
step is to prove that eq. ( 87 ) is true for the observable under consideration - with 



the same operator in the square brackets. 



5. Glasma fields, Long range rapidity correlations 
5.1. Glasma fields 



The factorization formula ( 90 1 indicates that in order to compute the single inclusive 
spectrum (or any other inclusive quantity, for which the same factorization holds) 
it is sufficient to compute the gluon spectrum in a classical field obtained by solving 
the Yang-Mills equation (72), and then to average over all the configurations of 
the color sources pi. 2- This classical field therefore plays a crucial role in the CGC 
description of heavy ion collisions, and it is worth spending some time discussing 
its properties. 

First of all, thanks to causality, one can divide space-time in four regions (see the 
left part of the figure 23). The region labeled zero is completely trivial: an observer 
in this region sees only the vacuum, and in this region the classical field is zero. 
Observers in the regions one or two see only one of the nuclei, but not the collision 



'This property is crucial : if this were not true, there would be a mixing of the sources of the two 
nuclei in the coefficients of the logarithms, preventing their factorization. 
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itself. In these regions, the solution of the classical Yang-Mills equations is also very 
simple. In light-cone gauge, it readsp^ 

A ± =0 , A\ 2 = -Ul 2 d i U 1 , 2 , (91) 
9 ' 

where U\ t 2 is a Wilson line constructed with the source p\ or pi respectively. In the 
regions one and two, the chromo-clcctric and chromo-magnetic fields are transverse 
to the collision axis. Note also that, since these fields are pure gauges, there is no 
field strength in the regions one and two. In other words, no energy is deposited in 
these region^ 

The most interesting region is the third one, because it causally connected to 
two nuclei. This region is the locus of all the events that happen after the collision. 
The gauge fields can be obtained analytically at r = + P^ but one has to resort 
to numerical methods beyond that. The problem is best formulated in the Fock- 
Schwinger gaug<£|±2iEl 

x + A~ + x~A + = , (92) 

and by exploiting the invariance under longitudinal boosts in collisions at very 
high energy. If we use the above gauge condition to parameterize the fields A ± as 
A^ 1 = ±x /?, then the function (3 and the transverse fields A 1 are independent of the 
rapidity r\. Solving numerically the Yang-Mills equations in the forward light-cone 
leads to the chromo-electric and chromo-magnetic fields shown in the right plot of 



the figure 23 At r = + , i.e. just after the collision, these fields have vanishing 



transverse components: the field lines are all parallel to the collision axis, and form 



•>This is of course expected, since all we have in these regions is a single (stable) nucleus moving 
at a constant speed. 

k In this gauge, the constraint of covariant current conservation [D^, = becomes trivial. 
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elongated tubular structures. Later on, as time increases, the transverse compo- 
nents of the fields become comparable to the longitudinal components. The g7asm^] 
designates these strong color fields that populate the system at early timesP^ 

In the longitudinal direction, the correlation length of the fields is infinite at 
leading order (since the system is boost invariant). When the leading log correc- 
tions are resummed, the distributions of sources evolve with rapidity, thus breaking 
the boost invariance. However, since the resummed terms are powers of a s A?7, it 
takes at least a rapidity shift A 77 ~ a J 1 in order to see an appreciable variation of 
the sources. Therefore, even after summing the leading logs, the glasma fields re- 
main coherent over rapidity intervals of order aj 1 . I n the transverse direction, the 
correlation length of the glasma fields is controlled by the saturation momentumf^Sl 
and is therefore of order Qj 1 . 

The fact that the field lines are parallel to the collision axis immediately after 
the collision leads to a very peculiar form for the energy-momentum tensor, 

T^ v = diag (e, e, e, — e) . (93) 

The most striking feature is that the longitudinal pressure is negative, which is 
reminiscent of strings stretching in the longitudinal direction. As we shall see in the 
section [6j this result is not the end of the story though, and further resummations 
are necessary due to the presence of instabilities in the solutions of the classical 
Yang-Mills equations. 



5.2. Long range rapidity correlations 

The correlation properties of the glasma fields have been proposed as the source 
of the long range rapidity correlations observed among pairs of hadrons in heavy 
ion collisions llSittUU The data is shown in the left part of the figure 24 It displays 



several features: a central peak (centered at A?7 = Aip = 0) that can be interpreted 
as due to collinear fragmentation of a fast particle, and a ridge-like structure, narrow 
in Aip and very elongated in A77. The rest of the discussion will focus on the latter. 

Before trying to interpret this correlation, one can reach a very general conclu- 
sion, based on causality, regarding the nature of the phenomena that may produce 
it. This is illustrated in the right part of the figure [24j Consider two particles A and 
B, and assume that they are correlated. They are detected long after the collision 
(compared to the strong interaction timescales), and their last interaction occurred 
on the freeze-out surface, at a proper time of order r/ « 10 fm/c. Between their 
last interaction and their detection, they traveled on straight lines, with an angle 
determined by their rapidity. From the point at which they crossed the freeze-out 
surface, draw a light-cone pointing in the past direction: these light-cones (in red 
and green in the figure) are the locus of the events that may have influenced A or 

'The word "glasma" is a contraction of glass (from colored glass condensate) and plasma (from 
quark-gluon plasma). 
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Au+Au central 




detection (—1 m/c) 

freeze out (-10 fm/c) 
latest correlation 



Fig. 24. Left: measured two-hadron correlation, as a function of the rapidity difference and 
azimuthal angle difference between the two hadronsPSH Right: causal structure of the two-particle 
correlations. 



B respectively. Any event outside these light-cones cannot possibly have any influ- 
ence, by causality. A correlation between A and B means that some event had an 
influence on both A and B; therefore this event must lie in the overlap of the two 
light-cones described above, i.e. in the region in blue in the figure. From this figure, 
we see that this overlap region extends only to a maximal time, 



olal.in 



< Tf e -|A»|/2 _ (94) 



This upper bound depends exponentially on the rapidity separation between the two 
particles, and therefore becomes very small for long range correlations in rapidity. 
This simple argument tells us that, regardless of their precise nature, the effects 
that produce these correlations must happen very shortly after the collisiorf"] or 
pre-exist in the wavefunctions of the colliding nuclei. 



The leading log factorization formula that generalizes eq. ( 90 ) to the two-particle 
correlation idHHHI 

(95) 



dN * fin n 1 w T 1 w T i dN ^ dNl 

[Dp, Dp 2 \ Wil/oJ W 2 [p, 



d 3 p 1 d 3 p 2 Leading Log J L 1 2J L lj L 2J d 3 p 1 d 3 p 2 L< 

Note that the integrand is the product of the single inclusive spectra for gluons of 



momenta p l and p 2 respectively, a consequence of eq. (85). This indicates that, at 
this order, all the correlations between the two particles result from the averaging 
over the sources pi,2, and therefore must come from the distributions Wifpi] and 
W 2 [p 2 ] themselves. 

A natural candidate for the formation of long range correlations in rapidity is the 
glasma color fields (112111121 since they are coherent over rapidity intervals that extend 
to at least A77 ~ aj 1 (in turn, the approximate boost invariance of the glasma fields 
is a consequence of the slow JIMWLK evolution of the distributions Wi,2, which is 



m This causality argument is very similar to the reason why, in cosmology, the observation of angular 
correlations in the Cosmic Microwave Background provides informations about the physics that 
prevailed long before the CMB was emitted. 
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an effect of order a s ). Note that the peculiar structure of the glasma field lines, that 
form longitudinal tubes at early times, is not essential to this argument. Another 
important parameter in the argument is the transverse size over which the glasma 
fields are coherent; this size is of order Qj 1 . Particles emitted from two distinct 
glasma flux tubes (i.e. separated by more than Qj 1 in the transverse direction) are 



not correlated (see the left figure 25 ) since they are produced by incoherent fields. 



The transverse size of the tubes controls the strength of the correlation. Indeed, two 





IF 



O 



Fig. 25. Left: particle emission from two distinct glasma flux tubes. Right: effect of radial flow 
on the angular distribution of the emitted particles. 



particles are correlated only if they come from the same tube. The probability for 
that is of order (Q S R)~ 2 , where R is the typical size of the transverse overlap in 
the collision (it coincides with the radius of the nuclei in a collision at zero impact 
parameter). 

However, the two-gluon correlation one gets from the glasma fields at early times 
has no angular dependence: the two particles are not correlated in azimuth. Note 
that there is no causality argument that tells that the azimuthal correlation must 
be created early; it can be produced later on by radial fJowIIilEIH The transverse 
pressure of the matter created in the collision sets it in motion radially, reaching 
speeds that are a sizable fraction of the speed of light. This collective radial motion 
transforms an azimuthally flat two-particle spectrum in a spectrum that has a peak 



around Aip — (see the right figure 24 ) . The width and height of this peak depend 
on the velocity of the radial flow. In the left plot of the fig ure [26 ] i s shown a crude 
calculation of the amplitude of this pealP^ (several work|ii£llHZI have performed 
more detailed studies of this effect), based on radial velocities extracted from the 
data itself.EH 

From the rapidity dependence of the distributions Wi[pi] and W 2 [p2], due to 
the JIMWLK evolution, it is possible to determine the dependence of the two- 



gluon spectrum on the rapidities of the two gluons, from eq. (95 1. Using some 



shortcuts (such as solving the mean field BK equation^] rather than the much more 

n The JIMWLK equation is considerably harder to solve numerically. The main idea is a reformu- 
lation of the equation as a Langevin equationpZHl that allows a numerical study of the JIMWLK 
equation on the lattice! 121 lThis study has shown that the Balitsky-Kovchegov equation is indeed a 
rather good approximation, at least as far as 2-point correlators are concerned. Recently, this direct 
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Fig. 26. Left: amplitude of the peak in Aip, compared to data. Right: estimate of the rapidity 
dependence of the correlation in the CGC (note: the central peak is due to j et fr agmentation and 
has been superimposed by hand) and comparison with data from PHOBOS.E21 

complicated JIMWLK equation) , one obtains the rapidity dependences^ shown in 
the right plot of the figure [26} 

6. Glasma evolution: isotropization and thermalization 
6.1. Resummation of the leading secular terms 

As we have seen in the previous section, the early stages of heavy ion collisions are 
dominated by strong color fields, the glasma. Initially, these chromo-electric and 
chromo-magnetic fields are purely longitudinal. A consequence of this is that the 
longitudinal pressure is negative (in fact, exactly the opposite of the energy density 
at r = + ). This is problematic in view of the many successes of the hydrodynamical 
description of the evolution of the matter produced in heavy ion collisions LL25HI32I 
because a large and negative pressure poses problems if used as initial condition for 
hydrodynamics . 

However, there is another fact that suggests that the CGC description of nucleus- 
nucleus collisions that we have described so far is still incomplete. It has been ob- 
served in several studies that classical solutions of the Yang-Mills equations suffer 
from instabilities J122HS21 that make them extremely sensitive to their initial condi- 



tion. The plot of the figure 27 shows the growth of such an unstable mode, when 
one disturbs a boost invariant classical solution by a small rapidity dependent per- 
turbation. These instabilities appear to be related to the well known Weibel in- 
stabilityp^HHSl or fomentation instability, in plasma physics. Many works have 
already investigated the possible role of these instabilities in the thermalization of 
the quark-gluon plasmaP^k^l 

In the CGC framework, these instabilities question the validity of the power 
counting that was the very basis for the organization of the expansion in power 
of g 2 . Indeed, this power counting implicitly assumes that if a classical field 

numerical approach has started to find its way into more phenomenological applic atinncl 122 l 12 3| 



42 Francois Gelis 



0.0001 

le-05 
le-06 

I 

L lc-07 
le-08 
. le-09 
| le-10 
le-11 
le-12 
le-13. 



— c +c, Exp(0.427 Sqrt(g fit)) 
c +e, Exp(0.00544 g 2 |i %) 




[Romatschkc, Venugopalan (2005)] 



1500 2000 2500 
g 1^* 



Fig. 27. Growth of some instabilities in classical Yang-Mills equations 

is of order g^ 1 and a small perturbation around it is initially of order 1, then 
the perturbation remains negligible compared to the background at all times. The 
existence of instabilities invalidates this assumption. It has been argued that the 
instabilities grow exponential!;^] as 

o( T ) ~ , (96) 

where the instability growth rate 7 is of the order of the saturation momentum. This 
means that it will become comparable in magnitude to the classical background field 
at a time 



7 



(97) 



is the time where the power counting rule in eq. ( 73 ) completely breaks down, 



and the one-loop corrections become as large as the leading order results. These 
terms, that are formally suppressed by powers of the coupling g, but with prefactors 
that grow with time, are called secular terms. 

In order to study the evolution of the system beyond this time, it is necessary to 
perform a resummation of the secular terms that have the fastest growth in time. A 



small amendment to the power counting rule of eq. ( 73 ) is sufficient to track these 
terms, 



loop - g 

The new rule is to assign a factor e 



T u ~ e^ T . (98) 
7T to each occurrence of the operator T. The 



reason for this new rule is eq. ( 86 1 , that shows that each power of T generates a 



perturbation on top of the classical background. 

In order to illustrate these rules, we have represented in the figure [28] some one- 
loop and two-loop contributions to the energy momentum tensor. With these rules, 
the one-loop graph on the left would bring a factor 

\ 2 



(ge^Y 



(99) 



°The square root inside the exponential is due to the longitudinal expansion of the system, that 
reduces the growth of the instability compared to a system that has a fixed volume. 
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Fig. 28. Examples of one-loop and two-loop contributions to the energy momentum tensor. The 
operators T are represented as blue dots. 



and the two-loop graph in the middle a factor 



/rn 



(100) 



while the two-loop graph on the right gives only 



(101) 



(These factors are all relative to the tree level result.) At the time r max , the first two 
contributions have the same order of magnitude as the leading order result, while 
the third contribution we have considered in this example is still suppressed. It is 
easy to resum all the terms that are leading at the time r max : this corresponds to 
all the graphs that have the same structure as the first two on the left of the figure 



28 with an arbitrary number of loops provided that these loops are not nested 



below the light-cone. The sum of these graphs can be obtained by exponentiating 



the operator in the square brackets in eq. (87). For a generic inclusive observable, 
this reads 



Ore 



exp 



G{u,v)T u T v 



(3(u)T % 



O, 



(102) 



By construction, this resummation contains the full LO and NLO results, and a 
subset of all the higher order terms, 



= o,n + a 



(103) 



full 



partially 



Eq. ( 102 ) can be turned in a much more practical expression. By noticing that the 



exponential of T is a translation operator, 



*^l-^ii 



itial 



exp 



initial I 



(104) 
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one can rewrite eq. ( |102[ ) as^] 



Crosum = J [DX] exp 



X(w) G 1 (u,v) x{v) 

,«6E 



LO ["4init+X + /3] • (105) 



In words, the resummed observable is obtained by shifting the initial value on E 
(in practice a surface of constant proper time tq) by a constant shift /3 and by a 
fluctuating shift \ that has a Gaussian distribution^] The variance of this Gaus- 
sian distribution can be calculated easily, and the functional integration over x 



in eq. ( 105 1 can be evaluated by a Monte-Carlo sampling. This resummation is 
sometimes described in the literature as classical statistical field theory. The same 
Gaussian average was also obtained in different approaches jiSiHUH A similar method 
has also been applied to cold atom physics, in problems related to Bose-Einstein 
condensationEHEZni 

In the rest of this section, we will present some numerical results illustrating the 
effect of this resummation. Because its implementation has not yet been completed 
for Yang-Mills theoryP^D we have tried it on a much simpler modeP^Hilil _ a rea [ 
scalar field with a (j> 4 coupling-, that shares with QCD some important features: 

• Scale invariance at the classical level in 3 + 1 dimensions. 

• Instabilities in the classical equations of motion, due to parametric reso- 
nancd 

Many aspects of this resummation can be studied by considering a system confined 
in a fixed volume, while other questions require to consider a system that expands 
longitudinally, as in a high energy collision. 

6.2. Evolution of a fixed volume system 

Consider first a system in a fixed volume. It is initialized with a large background 
field that mimics the glasma field of heavy ion collisions, to which we superimpose 



Gaussian fluctuations given by eq. (105). Then, each field configuration evolves 



according to the classical equation of motion, and observables are evaluated in its 



p In order to prove this formula, we need also 



r e -z z /2a 

f{x) = / dz f(x + z) 



q One can check that this Gaussian fluctuation of the initial field amounts to promoting the purely 
classical state >4i n it into a quantum mechanical coherent state centered at *4i n it- Coherent states are 
states that have the minimal extent in phase-space allowed by the uncertainty principle, and this 
extent is symmetrical between the coordinates and their conjugate momenta. Another equivalent 
point of view is that these fluctuations amount to filling each mode with 1/2 particle, thereby 

reproducing the ground state of a quantum oscillator. 

'Parametric resonance has been studied extensively in other situationsf5!5^Z!lE3 especially in 
inflationary cosmology. 
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solution. A Monte-Carlo sampling is performed in order to average over the Gaussian 
ensemble of fluctuations. 

Firstly, it is interesting to look at observables in a fixed loop order expansion, 
in order to highlight the pathologies caused by the presence of instabilities in the 



theory. This is illustrated in the figure 29 where we show the tree level (left) and 



one-loop (right) results for the energy density and the pressure. (The evolution 





Fig. 29. Energy density and pressure in <f> 4 scalar theory. Left: tree level result. Right: tree level 
+ one- loop result. 



really starts at t = 0. At t < 0, the system is initialized by a source that drives the 
background field to a large non-zero value.) One sees that the pressure oscillates 
in time - this means that there is no equation of state, since that would require 
a one-to-one correspondence between the pressure and the energy density. Even 
more problematic is the fact that at NLO, the oscillations of the pressure grow 
exponentially in time, making the predictions of the expansion in g 2 completely 
unreliable after a finite time. 

In contrast, the resummed pressure behaves in a very different way, as shown 
in the figure [30] The oscillations are damped, so that the pressure reaches a fixed 
value after a certain time. This asymptotic value of the pressure is related to the 
energy density by a very simple equation of state, 

e = 3P , (106) 

which is the expected relationship in a scale invariant theory in 3+1 dimensions. 
On similar timescales, the spectrum of excitations in the system also becomes con- 
siderably simpler, as shown by the spectral functions at two different times in the 
figure [3l| At the initial time, the spectral function has a complicated structure with 
more than one branch (including a branch that has a negative slope) . At later times 
(comparable to the relaxation time of the pressure), only one branch remains in the 
spectral function, suggesting that the system can be interpreted as a simple gas of 
quasiparticles. 

However, the relaxation of the pressure does not imply that the system has fully 
thermalized, as illustrated in the figure [321 Here, the background field is chosen such 
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that only a nonzero momentum mode has a large occupation number at the initial 
time, while all the other modes are empty. Very quickly, particles are produced in 
all the other modes, including the zero mode, mostly by elastic scatterings. At late 
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times, one can fit the occupation number by a function of the form, 

F ^ = ^~l (1 ° 7) 

where oj 2 , = k 2 + m 2 with m the quasiparticle mass obtained by studying the spec- 
tral function. T and /i are fitted to the numerical results. This functional form 
corresponds to the first two terms of the expansion of a Bose-Einstein distribution 
in the soft sectoiQ Interestingly, the best fit requires a nonzero chemical potential. 
Since the theory under consideration has no conserved particle number, this sug- 
gests that the system is not yet in chemical equilibrium even at the largest times 
considered. This is expected for the 4> A theory at weak coupling, because the inelas- 
tic cross-section (that can equilibrate the particle number) is much smaller than the 
elastic one. 

Another interesting feature is the fact that the fitted value of the chemical 
potential is very close to the mass of the quasiparticlcs. This, combined with the 
fact that there seem to be an excess of particles in the zero mode, suggests that the 
system undergoes Bose-Einstein condensation. We have confirmed this by showing 
that the occupation number in the zero mode is proportional to the volume of 
the system. Of course, in a theory where the particle number is not conserved, 
a Bose-Einstein condensate cannot be stable. Eventually, the inelastic collisions 
will eliminate the excess of particles that led to its formation, and dissolve the 
condensate. This is visible in the right plot of the figure [32] where the occupation 
number in the zero mode is shown as a function of time. After a very rapid growth, 
/(0) remains almost constant for a rather long time, to eventually decrease on 
even longer timescales (the timescales over which the inelastic processes become 
effective). As expected, by increasing the coupling, the condensate dissolves faster, 
because the inelastic collision rate is larger. 

Whether a Bose-Einstein condensate of gluons can form in heavy ion collisions 
is still an open question. The condition of overpopulation is initially also satisfied 
by the glasma fields b u t the inelastic rates are not parametrically suppressed 
compared to the elastic ones. Therefore it may happen that the condensate dissolves 
as fast as it forms, never really becoming a relevant feature of the time evolution. 



6.3. Evolution of a longitudinally expanding system 

In nucleus-nucleus collisions at high energy, another important question arises: is 
the longitudinal pressure comparable to the transverse pressure? This is an im- 
portant aspect in the applicability of hydrodynamics: if the two pressures are two 
different, then the viscous corrections are large. When discussing the glasma fields, 
we have seen that at leading order, the longitudinal pressure is in fact the opposite 



s The exponential tail of a Bose-Einstein distribution cannot be obtained in classical statistical 
field theory, because it is an approximation of the full theory which is valid only in the region 
where the occupation number is larger than unity. 
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of the energy density at r = + . This is in fact a generic feature in any system of 
longitudinally expanding fields. Indeed, energy and momentum conservation lead 
to the following equation (for a system homogeneous in the transverse plane) 

^ + ^=0. (108) 

OT T 

In order to a get a finite e at t = + , it is necessary that e + P L vanishes at r = + . 
As we have seen earlier, this is true in the glasma. This is also the case in the toy 
scalar model that we have considered in these numerical studies. 

A central question is whether the Gaussian fluctuations that are superimposed 
to the background field can cause the longitudinal pressure to become equal to the 
transverse one. This question can be addressed in the scalar </> 4 model considered in 
this section. We simply need to use r, rj coordinates, and to setup a boost invariant 
background field, i.e. independent of rj. The connection between this kind of lattice 




Fig. 33. Setup for the lattice study of a system in longitudinal expansion, as in a high energy 
collision. 



simulation and an actual collision is illustrated in the figure 33 It would be too 
costly to simulate the entire volume of the matter produced in a collision, therefore 
the lattice describes only a small patch of this matter that has a fixed extent in r\ 
and x± (and therefore stretches in z as time increases). 

The behavior of the energy density, and of the transverse and longitudinal pres- 
sures, is shown in the plot on the left of the figure [34] At very early times, the trace 
of the pressure tensor (2P T + P L ) has oscillations. These oscillations are quickly 
damped, and the system then obeys the equilibrium equation of state, 



e = 2P„ 



P, 



(109) 



However, at this stage of the time evolution, the transverse and longitudinal pres- 
sures still behave very differently. The longitudinal pressure decreases much faster, 
and even becomes negligible compared to the transverse pressure. This is explained 
naturally by the redshifting of the longitudinal momenta in an expanding system. 
This effect is illustrated in the right part of the figure [34] if one starts with a broad 
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Fig. 34. Left: time evolution of the energy density, and of the transverse and longitudinal pres- 
sures, in a longitudinal expanding system of scalar fields. Right: illustration of the redshifting of the 
longitudinal momenta due to longitudinal expansion. The black arrows represent the momentum 
vectors of particles, and the thin lines their trajectories if they move freely. 



distribution of longitudinal momenta in a given slice of rapidity 77, after some period 
of free streaming only the particles that have a momentum rapidity y equal to the 
space-time rapidity 77 remain in this slice. Thus, the longitudinal pressure (defined 
in the local rest frame of the matter, i.e. in a frame comoving along the rapidity 
slice) decreases rapidly. 

However, at some point, the longitudinal pressure stops decreasing. Instead, it 
increases exponentially and becomes very close to the transverse pressure, leading 
to an almost perfectly isotropic pressure tensor. This radical change of behavior of 
the longitudinal pressure seems related to the instabilities present in the system: at 
early times, the expansion rate (proportional to r _1 ) of the system is too large for 
the instabilities to be able to compete efficiently against the redshifting. It is only at 
later times, when the expansion rate has become low enough, that the instabilities 
become the driving force in the system. 

It is interesting to compare these results, obtained in classical statistical field the- 
ory, with hydrodynamical evolution. Naturally, since the transverse and longitudinal 
pressures are not equal, viscous corrections must be included in the hydrodynamical 
description to account for this difference. In a system which is homogeneous in the 
transverse plane and boost invariant, the simplest ansatz is to write 

p = - + ^ p =i-*L (no) 

3 3r ' L 3 3r ' [ ' 

where 77 is the shear viscosity of the system (this corresponds to first order hydro- 
dynamics). In an equilibrated scale invariant fluid, 77 is proportional to the entropy 
density s, and it is customary to quantify how viscous a fluid is by the dimension- 
less ratio rj/s. If we assume Stefan-Boltzmann's law for the entropy density, we have 
s s» e 3 / 4 , and since we know e, P T L as a function of r, we can extract from our results 
an effective value of 77/s. The result of this estimate is shown in the plot on the left 
of the figure 35 Clearly, the ratio 77/s is not constant. Its time dependence could be 
attributed to several causes: first order hydrodynamics (i.e. eq. (110)) may not be 
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Fig. 35. Left: estimate of the viscosity to entropy density ratio. Right: comparison between clas- 
sical statistical field theory and first order hydrodynamics, for various hydro initial times. 



a valid description of the system at these early times, and/or the system may still 
be too far from local equilibrium. Nevertheless, this estimate is instructive because 
it leads to an rj/s ratio which is much smaller than the perturbative estimate (blue 
band in the plot). This may be a manifestation of the anomalously small viscosity 
that has been conjecturedlHIIliS] for systems subject to instabilities and turbulence. 
Note also that the ratio we have extracted remains significantly above the value 
1/4-7T, obtained for certain gauge theories in the limit of infinite coupling. 

Another way to compare our results with first order hydrodynamics is to solve 
eq. (108) with the ansatz of eq. (110) for the longitudinal pressure, under the as- 
sumption that the ratio rj/s is a constant. To close the equation, we may assume 
again s ~ e 3 / 4 . A starting time To must be chosen for this comparison, and the differ- 
ence between P T and P L at this initial time determines the value of rj/s. For r > To, 
the evolution of the system is governed by eq. (108). This comparison is shown in 
the plot on the right of the figure 35 for several initial times. It appears that the re- 
laxation P L —> P T is much faster in the glasma description than in hydrodynamics, 
presumably due to the effect of instabilities. It is only if the hydrodynamical evolu- 
tion is started when the system is already well isotropized that it provides a good 
description of the evolution of the system. Naturally, eq. ( 108 ) is certainly valid 
because it is just a consequence of energy and momentum conservation. Therefore, 
this comparison questions the validity of the first order approximation to describe 
the viscous corrections (and the non-constancy of the ratio rj/s when extracted di- 
rectly from P T — P L points also in the same direction). The other conclusion one 
can make from this comparison, given the fact that hydrodynamics is the limit of 
kinetic theory when the mean free path goes to zero, is that instabilities are more 
efficient than collisions to isotropize the system. 



7. Conclusions 

In this review, we have presented the Color Glass Condensate effective theory, that 
describes the partonic content of a hadron or nucleus in the saturated regime. Gluon 
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saturation, that happens at high energy, is characterized by large gluon occupation 
numbers, and by non-linear effects not present in the dilute regime. In order to 
describe collisions involving hadrons in this regime, one needs a framework that 
provides a handle on multi-gluon Fock states in the hadron wavefunction, and the 
tools to perform calculations with these states. 

The Color Glass Condensate approximates the fast partons as stochastic classical 
color sources, while the usual field description is retained for the slow gluons. The 
requirement that observables be independent on the cutoff that separates these two 
types of degrees of freedom leads to a renormalization group equation (the JIMWLK 
equation) for the distribution of the classical sources. 

A crucial aspect of the physics of gluon saturation is that the non- linear dynamics 
generates a semi-hard scale, the saturation momentum. This scale, that increases 
with energy, plays a crucial role because it sets the value of the strong coupling 
constant, allowing a weak coupling treatment of gluon saturation. Nevertheless, 
even at weak coupling, the study of hadronic or nuclear collisions in the saturated 
regime is non-perturbative, because the smallness of the coupling is compensated 
by the large occupation numbers for the gluons below the saturation momentum. 

The Color Glass Condensate framework proves especially useful in the study of 
heavy ion collisions at high energy, because the bulk of particle production in these 
collisions is controlled by the physics of gluon saturation. Moreover, it has been 
shown that inclusive observables can be factorized (at leading logarithmic accuracy 
so far) as a convolution of universal distributions of color sources representing the 
gluon content of the two colliding nuclei, and an observable evaluated in the retarded 
classical field produced by these sources. 

The universality of these distributions of sources, and therefore of the classical 
color field they produce, confers to these field configurations -the glasma- a very 
central role in the CGC description of heavy ion collisions. At very short times after 
the collision, they are characterized by longitudinal chromo-electric and chromo- 
magnetic fields that are nearly boost invariant, with a transverse correlation length 
of the order of the inverse saturation momentum. In particular, these fields that 
have long range correlations in rapidity have been invoked in order to explain the 
observed two-hadron correlations in heavy ion collisions. 

The initial glasma fields form a system which is very from local equilibrium: it 
is a coherent state rather than a mixed thermal state, and it has a large negative 
longitudinal pressure. Moreover, these fields are unstable against small perturba- 
tions of their initial conditions. In fixed order CGC calculations, these unstable 
modes lead to secular divergences in quantities such as the energy-momentum ten- 
sor. This pathology can be cured by a resummation that amounts to letting the 
initial glasma field fluctuate around its classical value, with a Gaussian spectrum 
determined by a one-loop calculation. A numerical implementation of this resum- 
mation in a much simpler scalar field theory has recently shown that the instabilities 
lead to the isotropization of the transverse and longitudinal pressures shortly after 
the collision, and to full thermalization on longer timescales. Whether this is also 
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the case in QCD, and over which timescale, is the subject of ongoing works. 
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